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ABSTRACT 

In these companion papers, we study how the interrelated dynamics of sodium 
and potassium affect the excitability of neurons, the occurrence of seizures, and the 
stability of persistent states of activity. In this first paper, we construct a mathematical 
model consisting of a single conductance-based neuron together with intra- and 
extracellular ion concentration dynamics. We formulate a reduction of this model that 
permits a detailed bifurcation analysis, and show that the reduced model is a reasonable 
approximation of the full model. We find that competition between intrinsic neuronal 
currents, sodium-potassium pumps, glia, and diffusion can produce very slow and large- 
amplitude oscillations in ion concentrations similar to what is seen physiologically in 
seizures. Using the reduced model, we identify the dynamical mechanisms that give rise 
to these phenomena. These models reveal several experimentally testable predictions. 
Our work emphasizes the critical role of ion concentration homeostasis in the proper 
functioning of neurons, and points to important fundamental processes that may underlie 
pathological states such as epilepsy. 

Keywords: Potassium dynamics; bifurcation; glia; seizures; instabilities 
INTRODUCTION 

The Hodgkin-Huxley equations (Hodgkin and Huxley, 1952) have played a vital 
role in our theoretical understanding of various behaviors seen in neuronal studies both at 
the single-cell and network levels. However, use of these equations often assumes that 
the intra- and extra-cellular ion concentrations of sodium and potassium are constant. 
While this may be a reasonable assumption for the isolated squid giant axon, its validity 
in other cases, especially in the mammalian brain, is subject to debate. In this first of two 
companion papers, we investigate the role of local fluctuations in ion concentrations in 
modulating the behavior of a single neuron. 

Most studies investigating normal brain states have focused primarily on the 
intrinsic properties of neurons. Although some studies have examined the role that the 
extracellular micro-environment plays in pathological behavior (Bazhenov et al., 2004; 
Kager et al, 2000; Somjen, 2004; Park and Durand, 2006, Frohlich et al, 2008), little 
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attention has been paid to the cellular control of micro-environmental factors as a means 
to modulate the neuronal response (Park and Durand, 2006). 

In general, the intrinsic excitability of neuronal networks depends on the reversal 
potentials for various ion species. The reversal potentials in turn depend on the intra- and 
extracellular concentrations of the corresponding ions. During neuronal activity, the 
extracellular potassium and intracellular sodium concentrations ([K\ and [Na] h 
respectively) increase (Amzica et al., 2002; Heinemann et al, 1977; Moody et al, 1974; 
Ransom et al., 2000). Glia help to reestablish the normal ion concentrations, but require 
time to do so. Consequently, neuronal excitability is transiently modulated in a competing 
fashion: the local increase in [K\ Q raises the potassium reversal potential, increasing 
excitability, while the increase in [Na]i leads to a lower sodium reversal potential and 
thus less ability to drive sodium into the cell. The relatively small extracellular space and 
weak sodium conductances at normal resting potential can cause the transient changes in 
[K\ to have a greater effect over neuronal behavior than the changes in [Na]j, and the 
overall increase in excitability can cause spontaneous neuronal activity (Mcbain, 1994; 
Rutecki et al., 1985; Traynelis and Dingledine, 1988). 

In this paper, we examine the mechanisms by which the interrelated dynamics of 
sodium and potassium affect the excitability of neurons and the occurrence of seizure-like 
behavior. Since modest increases in [K] are known to produce more excitable neurons, 
we seek to understand ion concentration dynamics as a possible mechanism for giving 
rise to and perhaps governing seizure behavior. Using the major mechanisms responsible 
for the upkeep of the cellular micro-environment, i.e. pumps, diffusion, glial buffering, 
and channels, we mathematically model a conductance-based single neuron embedded 
within an extracellular space and glial compartments. We formulate a reduction of this 
model that permits a detailed analytical bifurcation analysis of the dynamics exhibited by 
this model, and show that the behavior produced by the reduced model is a reasonable 
approximation of the full model's dynamics. The effects of ion concentration dynamics 
on the behavior of networks of neurons is addressed in the companion article (Ullah et 
al., Submitted). 

Some related preliminary results have previously appeared in abstract form 
(Cressman et al, 2008). 
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METHODS 

1 . Full model 

Our full model consists of one single-compartment conductance-based neuron 
containing sodium, potassium, calcium-gated potassium, and leak currents, augmented 
with dynamic variables representing the intracellular sodium and extracellular potassium 
concentrations. These ion concentrations are affected by the neuron's intrinsic ionic 
currents as well as a sodium-potassium pump current, a glial current, and potassium 
diffusion. Finally, the concentrations are coupled to the membrane voltage equations via 
the Nernst reversal potentials. 

The conductance-based neuron is modeled as follows: 



dV 

dt 



c — = i Na +i K +i a 



I Na = ~g Na [™JV)fh(V-V Na )- g NaL (V - V Na ) 



(V-V K )- gKL (V-V K ) 



' n 4 , gAHpiCaV 
[* K \ + [Ca\ ) 

^Cl = ~SdL (V ~ Vci ) 

^ t =0[a g (V)(l-q)-/3 g (V)q], q = n,h 

^1 = -0.002g Ca (F - V Ca ) I {1 + exp(-(F + 25) / 2.5)} - [Ca\ 1 80 (1) 

The supporting rate equations are: 
mSV) = a m {V)l{a m {V)+P m {V)) 
a m {V) = 0. \{V + 30) /[l - exp(-0. l(V + 30))] 

/? m (D = 4exp(-(F + 55)/18) 

a n (V) = 0M(V + 34) /[l - exp(-0. \(V + 34))] 

J3„ (V) = 0. 125 exp( -(V + 44)/80) 

a h (V) = 0.07 exp(-(F + 44) /20) 

A(F) = l/[l + exp(-0.1(F + 14))] 

Note that the overall leak current consists of the final terms in the above 
expressions for I Na and I K , plus I C1 ; similar leak currents were used by (Kager et al., 
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2000). Also, the gating variable m is assumed to be fast compared to the voltage change; 
we therefore assume it reaches its equilibrium value m x immediately (Rinzel, 1985; 

Pinsky and Rinzel, 1994). Finally, the active internal calcium concentration is used only 
in conjunction with the calcium-gated potassium current in order to model the adaptation 
seen in many excitatory cells (Mason and Larkman, 1990; Wang, 1998). 

The meaning and values of the parameters and variables used in this paper are 
given in Table 1 . 

The potassium concentration in the interstitial volume surrounding each cell was 
continuously updated based on K + currents across the neuronal membrane, Na + -K + 
pumps, uptake by the glial network surrounding the neurons, and lateral diffusion of K + 
within the extracellular space. Thus, we have 

^ = -0.337, -Ifil^ -I gUa -V (2) 

The factor 0.33mM.cm /ucoul converts current density to rate-of-change of concentration 
(see Appendix A). The factor f3 corrects for the volume fraction between the interior of 
the cell and the extracellular space when calculating the concentration change and is 
based on Mazel et al. (1998), McBain et al. (1990), and Somjen (2004). 

The sodium-potassium pump is modeled as a product of sigmoidal functions as 
follows: 



^ pump 



1.0 



l.0 + ex V ((25.0-[Nal)/3.0) 



1.0 + exp(5.5-[*] o ), 

Normal resting conditions are attained when p= 1.25mM/sec. Each term saturates 
for high values of internal sodium and external potassium, respectively. More 
biophysically realistic models of pumps, such as those in (Lauger, 1991) produce 
substantially similar results. 

The capacity of glial cells to remove excess potassium from the extracellular 
space is modeled by 

G , 



1.0 + exp((18-[if| o )/2.5) 

This highly simplified model incorporates both passive and active uptake into a single 
sigmoidal response function that depends on the extracellular potassium concentration 
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alone. Normal conditions correspond to G g n a = 66mM/sec, and [K] = 4.0mM. A similar 
but more biophysical approach was used in (Kager et al, 2000). Two factors allow the 
glia to provide a nearly insatiable buffer for the extracellular space. The first is the very 
large size of the glial network. Second, the glial endfeet surround the pericapillary space, 
which, through interactions with arteriole walls, can effect blood flow; this, in turn, can 
increase the buffering capability of the glia (Paulson and Newman, 1987, Kuschinsky et 
al, 1972, McCulloch et al, 1982). 

The diffusion of potassium away from the local extracellular micro-environment 
is modeled by 

diff VL J o o,x' 

Here, k g x is the concentration of potassium in the largest nearby reservoir. 

Physiologically, this would correspond to either the bath solution in a slice preparation, 
or the vasculature in the intact brain (noting that [K] G is kept below the plasma level by 
trans-endothelial transport). For normal conditions, we use k o oo = 4.0 mM. The diffusion 

constant s , obtained from Fick's law, is s=2D/ Ax 2 , where we use D = 2.5x10" 
cm /sec for K in water (Fisher et al., 1976) and estimate Ax « 20 /um for intact brain 
reflecting the average distance between capillaries (Scharrer, 1944); thus s= 1.2Hz. 

To complete the description of the potassium concentration dynamics, we make 
the assumption that the flow of Na + into the cell is compensated by flow of K + out of the 
cell. Then [K\j can be approximated by 

[K\ =140.0mM+(18.0mM-[7Yfl] ; ), (3) where 

140.0 mM and 18.0 mM reflect the normal resting [K\i and [Na] , respectively. The 
limitations of this approximation will be addressed in the discussion section. 

The intra- and extracellular sodium concentration dynamics are modeled by 

d[Na\ L, 

— — ^ = 0.33^-3/ (4) 

dt P pump 

[Na\ a = l44.0mM-j3([Nal-lS.0mM). (5) 
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In equation (5), we assume that the total amount of sodium is conserved, and hence only 
one differential equation for sodium is needed. Here, 144. OmM is the sodium 
concentration outside the cell under normal resting conditions for a mammalian neuron. 

Finally, the reversal potentials appearing in equation (1) are obtained from the ion 
concentrations via the Nernst equation 



V m = 26.64 In 



[Na\ 



V K =26.64 In 



V a =26.64 In 



[*1 . 

m 



With the leak conductances listed above, the chloride concentrations were fixed at 
[Cl\ =6.0mM and [C/] o =130mM. 

Thus, the dynamic variables of the full model are V, n, h, [Ca] , [K] o , and [Na\ . 

Despite the fact that we have neglected many features of real mammalian cells (such as 
the geometrically complex dendritic and axonal structure and the related spatially 
complex distribution of channels and cotransporters, as well as the presence of immobile 
anions which are strictly required to maintain electric and osmotic balance), our model 
captures the essential dynamics that we wish to explore. 

In the results section, we will be interested in varying the parameters G lia , s , 

k o x , and p . We will present our results in terms of parameters normalized by their 

normal values, for example, G gUa — G glia I G glia normal , where the overbar indicates the 
normalized parameter. 



2. Reduced model 

In order to more effectively study the bifurcation structure of the model presented 
above, we formulated a reduction by eliminating the fast-time-scale spiking behavior in 
favor of the slower ion concentration dynamics. This is accomplished by replacing the 
entire Hodgkin-Huxley mechanism with empirical fits to time-averaged ion currents. 
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Using the membrane conductances from the full model, we fixed the internal and external 
sodium and potassium concentration ratios and allowed the model cell to attain its 
asymptotic dynamical state, which was either a resting state or a spiking state. Then, the 
sodium and potassium membrane currents were time-averaged over one second. These 
data were fit to products of sigmoidal functions of the sodium and potassium 
concentration ratios, resulting in the (infinite-time) functions 

I Na ^Na\ l[Na\,[K\ /[K\) and I K ^Na\ I [Na\,[K\ I [K]). Details are available 

in Appendix A. I Nax is nearly identical to I Kx , differing significantly only near normal 
resting concentration ratios due to differences in the sodium and potassium leak currents. 

Thus, our reduced model consists of equations (2-5), with I Na and I K replaced 
with the empirical fits described above (see Appendix A for additional details). 

3. Bifurcations in the reduced model 

Our main results in this paper consist of identifying bifurcations in the reduced 
model and analyzing their implications for the behavior of the full model. We observe 
that depending on the various parameters, the ion concentrations in the reduced model 
approach either stable equilibria, and thus remain constant, or they approach stable 
periodic orbits, and thus exhibit oscillatory behavior (Fig 1). As parameters are changed, 
the stability of these solutions change. This happens through bifurcations (a good general 
reference is Strogatz (1994)). Most relevant for our purposes are the Hopf bifurcation, 
and the saddle-node bifurcation of periodic orbits. In a Hopf bifurcation, an equilibrium 
solution either gains or loses stability, and simultaneously, a periodic orbit either appears 
or disappears from the same point 1 . In a saddle-node bifurcation of periodic orbits, a pair 
of periodic orbits - one stable and one unstable - either appears or disappears suddenly, 
as if out of "thin air". 

Bifurcation diagrams were obtained using XPPAUT (Ermentrout, 2002). Code for 
both of our models is available from ModelDB. 



1 Depending on the stabilities of the equilibrium and the periodic orbit involved, Hopf bifurcations are 
classified as sub- or supercritical. 

2 http://senselab.med.yale.edu/modeldb/ 
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RESULTS 



1. Overview 

In an experimental slice preparation, an easily-performed experimental 
manipulation is to change the potassium concentration in the bathing solution. Such 
preparations have been used to study epilepsy (Jensen and Yaari, 1997; Traynelis and 
Dingledine, 1988; Gluckman, et al. 2001). At normal concentrations (~4mM), normal 
resting potential is maintained. However, at higher concentrations (8mM, for example) 
bursts and seizure-like events occur spontaneously. 

We begin discussing the dynamics of our models by considering a similar 
manipulation, corresponding to varying the normalized parameter k g x . In the full model, 

setting k = 2.0 (i.e., doubling the normal concentration of potassium in the bath 

solution) leads to spontaneously-occurring prolonged periods of rapid firing, as illustrated 
in the top trace of Fig. 1. These oscillations are remarkably similar to experimental results 
reported by several investigators (see, for example, Figures 1 and 6 of Jensen and Yaari 
(1997), in which the authors use a high potassium in vitro preparation, and Figure 2 of 
Ziburkus et al. (2006), which reports results from a 4-aminopyridine in vitro preparation). 
Each event lasts on the order of tens of seconds and consists of many spikes, each of 
which occurs on the order of 1 ms. Thus, the full model contains dynamics on at least two 
distinct time scales that are separated by four orders of magnitude: fast spiking from the 
Hodgkin-Huxley mechanism, and a slow overall modulation. The solid traces in the 
middle and bottom panels show that this slow modulation corresponds to slow periodic 
behavior in the sodium and potassium ion concentrations, respectively. 

Our reduced model was constructed in order to remove the fast Hodgkin-Huxley 
spiking mechanism and focus attention on the slow dynamics of the ion concentrations. 
The dashed traces in the middle and bottom panels of Fig. 1 show the sodium and 
potassium ion concentrations obtained from the reduced model for the same parameters 
used above. Although these traces are not identical to those of the full model, it is evident 
that the reduced model captures the qualitative behavior of the ion concentrations quite 
well. 
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The separation of time scales achieved by our model reduction (see, for example, 
Rinzel and Ermentrout (1989); Kepler, et al, (1992)) yields a model that is amenable to 
numerical bifurcation analysis. Knowledge of the bifurcation structure in turn informs us 
about the dynamical mechanisms that underlie the full model. In the following, we will 
first describe the main dynamical features of the reduced model, and then examine the 
implications for the behavior of the full model. 

2. Analysis of the reduced model 

Fig. 2 shows a bifurcation diagram obtained using the reduced model. This 
diagram plots the minimum and maximum asymptotic values of the extracellular 
potassium concentration [K] o versus a range of values of the reservoir's normalized 

potassium concentration k g x . For low values of this parameter, [K] o is observed to settle 
at a stable equilibrium. The value of [K] o corresponding to this equilibrium increases 
with k g x until the equilibrium loses stability via a subcritical Hopf bifurcation at 
k ox « 1.9 . This means that at this point, an unstable periodic orbit collapses onto the 
equilibrium, and both disappear. [K] o is subsequently attracted to a coexisting large- 
amplitude stable periodic orbit . Thus, large-amplitude oscillations in [K] o appear 
abruptly. These oscillations persist as k g x is increased until the same sequence of 
bifurcations occurs in the opposite order at k « 2. 13 . At this higher value of k , the 

1 1 0,00 0,00 

unstable equilibrium undergoes a subcritical Hopf bifurcation, becoming stable and 
giving rise to an unstable periodic orbit whose amplitude quickly rises with increasing 
k g x . This orbit then collides with the large-amplitude stable periodic orbit at k x « 2.15, 

and both orbits disappear in a saddle -node bifurcation of periodic orbits. In this manner, 



3 The stable and unstable periodic orbits involved in this scenario appear via a saddle-node bifurcation at a 
slightly smaller parameter value that is extremely close to that of the Hopf bifurcation. Thus, the sequence 
of bifurcations is not immediately apparent in Fig. 2. The abruptness of these transitions, and the difficulty 
in resolving them numerically, is due to the "canard" mechanism (Dumortier and Roussarie, 1996; 
Wechselberger, (2007)). 
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the periodic behavior of [K] o is terminated. For still higher values of k o x , [K] o 

approaches the equilibrium values shown at the far right in Fig. 2. 

In order to examine the boundaries of the oscillatory behavior described above 
with respect to pump strength p , the diffusion coefficient s , and glial buffering strength 

G glia , we constructed the bifurcation diagrams shown in Fig. 3. First, we fixed all 

parameters at their normal resting values except k o x , which we set to 2 in order to obtain 

the oscillatory behavior discussed above. We then separately varied p , s , and G glia 
away from their normal resting values. If p and s are increased from their nominal 
values of 1 (Fig. 3a, b), we see that the oscillatory behavior terminates in a manner 
similar to that described above; that is, an unstable periodic orbit appears via a subcritical 
Hopf bifurcation which grows until it collides with and annihilates the stable periodic 
solution at a saddle-node bifurcation of periodic orbits. (This is most apparent in Fig 
3(a).) The same scenario applies as these parameters are decreased 4 , 5 . The situation is 
similar for the glial strength G gUa , except that on the left, we see no saddle-node 

bifurcation of periodic orbits for positive values of G glia (Fig. 3 c). 

It is notable that if s or G glia are reduced from their normal values, larger 

amplitude oscillations in [K] o occur. This is because, in both cases, the trafficking of 

potassium away from the extracellular space is impeded, and consequently, [K] o builds 

up more effectively during the spiking phase of the cell's activity (see Fig. 6, below). In 
contrast, changing p in either direction results in very little change in the amplitude of 

the [K] o oscillations. Furthermore, the bistable region on the right side of the p diagram 
is quite wide, and hence hysteretic behavior as p is varied across this region may be 
particularly amenable to experimental observation (e.g. with ouabain). Finally, we note 
that the right branch of stable equilibria in Fig. 3(a) corresponds to higher values of 
[K] o than the left branch. This is because the cell is active - either spiking or in 

4 A canard similar to that described previously occurs here, so that the Hopf and the saddle -node 
bifurcations on the left sides of Figs. 3 (a) and (b) occur in extremely narrow intervals of the parameter. 

5 In Fig. 3 (a), the equilibrium curve does not extend all the way to zero because of the constant chloride 
leak current. 
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depolarization block - and thus there is a large membrane potassium current flowing into 
the extracellular space which must be balanced by pumps and otherwise "normal" 
diffusion and glial currents. 

The two-parameter bifurcation diagram shown in Fig. 4 provides a more complete 
understanding of the oscillatory behavior of [K] o in our reduced model with respect to 

the variation of both £~and G .. , with k = 2. The dashed lines at G ,,. = 1 and s = \ 

gua 7 o^cc glia 

correspond to the one-dimensional bifurcation diagrams shown in Figs. 3b and c. The 
solid curves represent Hopf bifurcations, and the intersection of these dashed lines with 
the Hopf curves correspond to the Hopf bifurcations (points) in the earlier figures. Thus, 
the Hopf curves define a region within which [K] o is obliged to oscillate, because the 

only stable attractor is a periodic orbit. To facilitate discussion, we refer to this region as 
the "region of oscillation", or RO. Outside of the RO, stable equilibrium solutions for 
[K\ exist 6 . 

The dashed line at G glia = 1 .75 corresponds to the one-dimensional bifurcation 

diagram shown in Fig. 5. This is drawn at the same scale as Fig. 3b (the G glia = 1 

diagram) to facilitate comparison. We note that the amplitude of the oscillation in [K\ is 
significantly smaller in this region of the RO. Furthermore, the Hopf bifurcation on the 
right (at about s= 3.2) is now supercritical. This means that the amplitude of the [K\ 
oscillation decays smoothly to zero as this point is approached from the left. 

3. Analysis of the full model 

We now investigate whether the dynamical features identified above in our 
reduced model correspond to similar features in our full model. Fig. 6 shows traces of the 
membrane voltage (upper traces), [K\ (solid lower traces), and [Na\ (dashed lower 

traces) versus time, all obtained from the full model, corresponding to the parameter 
values marked by the numbered squares in Fig. 4. For the regions outside of the RO, the 
reduced model predicts stable equilibrium solutions for the ion concentrations. For 
example, at point 1, [K\ is slightly elevated at a value near 6mM, and the membrane 

6 Note that oscillations may persist slightly outside of the RO, where a stable periodic orbit coexists with 
the stable equilibrium solution; see, for example, the right side of Fig. 3(a). 
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voltage of the full model remains constant at -62 mV (not shown). However, at point 2, 
the full model exhibits tonic firing, as shown in Fig. 6a. Here, [K\ is sufficiently high 
such that the neuron is depolarized beyond its firing threshold, and [K\ remains 
essentially constant with only small perturbations of order 0.1 mM due to individual 
spikes (Fig. 6a, lower panel) (Frankenhaeuser and Hodgkin, 1956). (These spike 
perturbations disappear in the reduced model, as they are smoothed-out by the averaging 
in our model reduction procedure.) 

Within the RO, the reduced model predicts periodic behavior with relatively large 
and slow oscillations in the ion concentrations. Points 3-7 in Fig. 4 correspond to Figs. 6 
(b-f), which show various bursting behaviors of the full model. To facilitate comparison, 
the time scale for these figures (Fig. 6 b-f) is the same, showing 100 seconds of data. In 
addition, the voltage and concentration scales are also the same, except for the 
concentration scale in (b). We make the following observations from Fig. 6. As the RO is 
traversed from low to high a (keeping G glia fixed) in Fig. 4, the bursts become more 

frequent (compare Fig. 6c, d; points 4, 5 and Fig.6e, f; points 6, 7). In addition, the shape 
of the burst envelope changes due to the decreasing amplitude of the [K\ oscillations. 
Note also in Fig. 6b (point 3) that the peak of the [K\ concentration is nearly 40mM, 
large enough to cause the neuron to briefly enter a state of quiescence known as 
"depolarization block" (see companion paper, Ullah, et al., Submitted). Finally, the 
within-burst spike frequency is essentially constant in Figs. 6 (b)-(f); it does peak in 
concert with the peak of [K] o , however, this effect is weak. 

Returning briefly to the reduced model, we address in Fig. 7 how the location of 
the RO changes with k as all other parameters are kept at their normal values. As k 

G o,oo 1 1 0,00 

is increased from its normal value at 1 , the RO emerges around a value of k - 1 .77 , as 

represented by the gray line on the left of Fig. 7a. (This observation is consistent with 
Rutecki , et al. (1985)). As k g x increases, the right edge of the RO shifts towards the 

right, crossing the normal values of ^Tand G glia at (1,1) when k oa = 1.9, as shown by the 

thick solid line in Fig. 7a. As k is further increased to 2.0, normal conditions for glial 

pumping and diffusion are well inside the RO as shown in Fig. 7b, and correspondingly, 
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we observe bursting/seizing behavior in the full model. Fig. 7c shows the RO when the 
reservoir concentration has been increased to 2.1. Here the left side edge of the RO is 
just about to cross the point (1,1), after which the ion concentrations assume stable 
equilibrium values in the reduced model. Note that these stable equilibria actually 
correspond to a state of rapid tonic firing in the full model, as in Fig. 6(a); for much 
higher values of k g x (greater than approximately 7.2mM), [K\ remains constant in the 

reduced model, but the full model eventually gives rise to depolarization block. 

We conclude by reporting other kinds of bursting behaviors seen in our full 

model. Fig. 8a shows a time trace of the membrane voltage for k -6, G , -OA and 

s = 0.4 . These bursts are fundamentally different from those shown in Fig. 6. In 
particular, the extracellular potassium concentration is quite elevated, and thus the 
periods of quiescence correspond not to resting behavior, but rather to a state of 
depolarization block. In addition, the bursts themselves have a rounded envelope, as 
opposed to the (approximately) square envelopes of the events shown in Fig. 6. This 
behavior is consistent with the experimental observations of (Ziburkus et al., 2006), in 
which interneurons were seen to enter depolarization block and thus give way to 
pyramidal cell bursts. Bikson et al. also observed depolarization block in pyramidal cells 
during electrographic seizures (see Figure Id from Bikson et al, 2003). We have also 
experimentally observed (in oriens interneurons exposed to 4-aminopyridine) relatively 
continuous "burst" firing without any quiescent intervals, as seen in Fig.8b. In this figure, 
the neuron fires continuously, but with a wavy envelope due to the oscillating ion 
concentrations. We include these patterns to complete the description of the repertoire of 
single cell bursting behaviors seen in our models. 

DISCUSSION 

We have created a model which can be used to investigate the role of ion 
concentration dynamics in neuronal function, as well as a reduced model which is 
amenable to bifurcation analysis. Such bifurcations indicate major qualitative changes in 
system behavior, which are in many ways more predictive and informative than purely 
quantitative measurements. In particular, we show that under otherwise normal 
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conditions, there exists a broad range of bath potassium concentration values which 
yields seizure-like behavior in a single neuron that is both qualitatively as well as 
quantitatively similar to what is seen in experimental models (Ziburkus et al., 2006; Feng 
and Durand, 2006; McBain, 1994). In fact, the values of extracellular concentration used 
in those experiments are quantitatively consistent with the range of concentrations shown 
here to exhibit seizures. Furthermore, the stable periodic oscillations in the extracellular 
potassium concentration which result from varying various experimentally and 
biophysically relevant parameters suggests that these effects may be an important 
mechanism underlying epileptic seizures. 

In formulating our models, we made several approximations. The two most 
severe, insatiable glial buffering and the assumption that internal potassium can be 
calculated using equation (3), should hold for times that are long compared to the time 
scale of individual spikes, bursts, and seizures. However, for even longer times (on the 
order of thousands of seconds), the saturation of glia as well as the decoupling of the 
internal potassium from the sodium dynamics will lead to more substantial errors in the 
calculated results. It is possible that glia do not fully saturate, if, as suggested by Paulson 
and Newman et al (1987), the glia siphon excess potassium into pericapillary spaces via 
the astrocyte network. Nevertheless, one can understand a slow partial saturation of the 
glial network as a slow decrease in G glia , for example. Consequently, the system may 

enter or leave parameter-space regions in which oscillating ion concentrations exist (e.g., 
see Figs. 3(c) and 4). This long-term behavior could be used to more accurately model 
the temporal dynamics of the glial siphoning system. 

It is important to note the limitations of our models with respect to extremely long 
time scales. If the reservoir concentration k o x remains only slightly elevated for a long 

period of time the model cell will ultimately reach a new fixed point in [K\ nearly equal 
to the bath concentration. A stable resting state should exhibit some degree of robustness 
in its micro-environment. However, as the system drifts further from the normal state, we 
should not expect such homeostasis to persist; the internal potassium will in general drift 
to higher or lower values depending on the wide variety of pumps, cotransporters, or 
channels inherent to the cell. When the internal potassium is integrated separately (not 
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shown here), we see both upward and downward drift depending on model parameters as 
well as the initial conditions for the ionic concentrations. Therefore in our model the 
seizure-like events, as well as the tonic firing reported here, appear to be transients on 
extremely long time scales. Of course, such phenomena are also transients on long time 
scales in animals and humans. 

Although our reduced model does a good job reproducing the qualitative results 
of our full spiking model, there are regions were the two models disagree. These two 
models produce very good agreement in the region of the two-parameter graph presented 

in Figure 4. However the reduced model predicts the cessation of all oscillations as G gUa 

is increased past a value of 4 (not shown), whereas the full model exhibits burst-like 
oscillations for far greater values. This discrepancy is due to the use of relatively simple 
functions used in our fitting of the time-averaged Hodgkin-Huxley currents (see 
Appendix A). A more sophisticated fit of these data would improve the agreement 
between our two models. 

Our work points out the important role that ion concentration dynamics may play 
in understanding neuronal dynamics, including pathological dynamics such as seizures. 
Of course, in realistic situations, these are due to a combination of local environmental 
conditions and electrical and chemical communication between cells (see the 
accompanying paper, Ullah, et al, Submitted). The models presented here, however, 
demonstrate that recurring seizure-like events can occur in a single cell that is subject to 
intra- and extracellular ion concentration dynamics (see also discussion in (Kager, et al., 
2000, Kager, et al., 2007) regarding single-cell seizure dynamics). In addition, we have 
identified the basic mechanism that can give rise to such events: Hopf bifurcations that 
lead to slow oscillations in the ion concentrations. 
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FIGURE LEGENDS 

Fig. 1: Comparison of the reduced model to the full spiking model. The top plot shows 
the membrane voltage for the neuron in the full model. The middle traces show [K] o for 
both the full model (solid line) and the reduced model (dashed line). The bottom traces 
show [A/a], with the same convention. All data were integrated with an elevated bath 

potassium concentration at k ox — 2.0 , with all other parameters set to their normal 
values. 



Fig. 2: The bifurcation diagram for [K] o as a function of the bath potassium 

concentration k a x , revealing a region of oscillatory behavior. All other parameters were 

set equal to their normal values. Triangles represent equilibria (i.e., steady states), and 
circles represent periodic orbits (i.e., oscillatory behavior); stability is denoted by solid 
(stable) and open (unstable) symbols. 

Fig. 3: Bifurcation diagrams for [K] o as a function of (a) the normalized pump strength, 
(b) the diffusion coefficient, and (c) the glial strength. All plots were produced with an 
elevated bath potassium concentration at k — 2.0 . 

1 n.cc 



Fig. 4: A two-dimensional bifurcation diagram showing the boundaries of the region of 
oscillation (RO) as a function of the diffusion coefficient and the glial strength. The black 
curves denote Hopf bifurcations; within this region, the ion concentrations exhibit 
oscillatory behavior. The dashed lines correspond to the one-dimensional bifurcation 
diagrams in Figures 3b, 3 c and 5 (see text). Examples of the dynamics of the full model, 
obtained at parameter values corresponding to the numbered points, are shown in Figure 
6. 

Fig. 5: The one-dimensional bifurcation diagram for [K] o as a function of the normalized 
diffusion coefficient a for£ =2.0 and G ,,. = 1.75. 

0,00 glia 
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Fig. 6: Examples of the dynamics of the full model, obtained at parameter values 
corresponding to the numbered points in Figure 4. The top trace shows the membrane 
voltage and the lower traces show [K\ (solid trace) and [Na]i (dashed trace) on the same 
time scale. 

Fig. 7: The effect of changing the bath concentration on the location of the region of 
oscillation (RO) is illustrated for p = 1 and various values of the bath potassium 

concentration k a x . The square represents normal values of the diffusion and glial 
strength. In (a), the RO is seen to appear and move to the right as the bath potassium 
concentration is increased from k g x - 1 .77 (grey curve) to k =1.9 (black curve), 

where it intersects the square (compare the left bifurcation in Figure 2). In (b), the normal 
square lies within the RO for k = 2.0 . In (c), k = 2. 1 , the RO has moved further to 

A o,oo v /J 0,00 J 

the right, and the square is close to the left boundary. 

Fig. 8: Other bursting patterns. Traces similar to those in Figure 6 obtained with the full 
model for (a) k = 6.0, G = 0.1, and s = 0.4 ; and (b) same, but with k reduced 

v 7 o,oo 7 gtia \ / Q^oo 

slightly. The quiescent states in (a) correspond to depolarization block; see text for 
further description. 
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Table 1 : Model variables and parameters 



Variable 



Units 



Description 



V 
Im 
h 
h 

mJ[V) 

h 

n 

a(V) 

P(V) 

[Cal 

V Na 

V K 

[Na] 
[Nal 
[K] 
[K\t 

Ipump 



mV 

uA/cm 

uA/cm 

2 

uA/cm 



mM 
mV 

mV 

mM 

mM 

mM 

mM 

mM/sec 

mM/sec 

mM/sec 



Membrane potential 

Sodium current 

Potassium current 

Leak current 

Activating sodium gate 

Inactivating sodium gate 

Activating potassium gate 

Forward rate constant for transition between 

open and close state of a gate 

Backward rate constant for transition between 

open and close state of a gate 

Intracellular calcium concentration 

Reversal potential of persistent sodium 

current 

Reversal potential of potassium current 
Extracellular sodium concentration 
Intracellular sodium concentration 
Extracellular potassium concentration 
Intracellular potassium concentration 
Pump current 

Potassium diffusion to the nearby reservoir 
Glial uptake 



Parameter 



Value 



Description 



C 

gNa 

gK 

gAHP 

gKL 

gNaL 

gCIL 

Va 

gCa 
V Ca 

P 
P 

GgUa 
S 

[Cl\i 
[CQq 



l|uF/cm 

100mS/m 2 

40mS/m 2 

0.01mS/m 2 

0.05mS/m 2 

0.0175mS/m 2 

0.05mS/m 2 

3sec 

-81.93mV 
O.lmS/m 2 
120mV 
7.0 



1.25mM/sec 

66mM/sec 

1.2sec _1 

4.0mM 

6.0mM 

130.0mM 



Membrane capacitance 

Conductance of persistent sodium current 

Conductance of potassium current 

Conductance of afterhyperpolarization current 

Conductance of potassium leak current 

Conductance of sodium leak current 

Conductance of chloride leak current 

Time constant of gating variables 

Reversal potential of chloride current 

Calcium conductance 

Reversal potential of calcium 

Ratio of intracellular to extracellular volume 

of the cell 

Pump strength 

Strength of glial uptake 

Diffusion constant 

Steady state extracellular concentration 
Intracellular chloride concentration 
Extracellular chloride concentration 
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APPENDIX A 

Current to concentration conversion: 

In order to derive the ion concentration dynamics, we begin with the assumption 
that the ratio of the intracellular volume to the extracellular volume is ft = 7.0 (Somjen, 
2004). This corresponds to a cell with intracellular and extracellular space of 87.5% and 
12.5% of the total volume respectively. For the currents across the membrane, 
conservation of ions requires 

Ac Vol. = -Ac Vol , 

11 o o 

where c and Vol represent ion concentration and volume respectively, A indicates change, 
and the subscripts i, o correspond to the intra- and extracellular volumes. The above 
equation leads to 

( Vol } Ac 

Ac. = -Ac 



Vol J J3 

Let / be the current density in units of uA/cm from the Hodgkin-Huxley model. 
Then, the total current i = IA entering the intracellular volume produces a flow of 

charge equal to AQ = i total At in a time At , where A is membrane area. The number of ions 
entering the volume in this time is therefore AN = i total At I q where q is 1 .6xl0" 19 coul. 
The change in concentration Ac. = AN I N A Vol. depends on the volume of the region to 
which the ions flow, where Avogadro's number Na converts the concentration to molars. 
The rate of change of concentration, or concentration current dc. I dt-i . , is related to 

° ! C,l ' 

the ratio of the surface area of the cell to the volume of the cell as follows 

Ac. AN i ft , IA I 

j i_ total 

e '' ~ At ~ AtVol.N A ~ qVol.N A ~ qVol.N A ~ a ' 

2 

For a sphere of radius 7um, a = 21mcoul/M.cm . An increase in cell volume 
would result in a smaller time constant and therefore slower dynamics. 

For the outward current the external ion concentration is therefore given as 

i =Bi =^- = 0.33/ 

C,0 " C.l 

a 
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Equations for reduced model: 

The reduced model uses empirical fits of the average membrane currents of the Hodgkin- 
Huxley model neuron, as described in the main text. The fits are given below. 

K o/l =[K] o /[Kl 
Na Uo =[Na\l[Nd\ o 



g x = 420.0(l-A 1 (l-B 1 exp(- / / 1 Na i/o )) 1 ") 
g 2 = explo-^l.O-^K^Kl.O+exp^Na^))) 
g 3 = (l/(l+exp(o- 3 (1.0+// 3 Na i/o a 3 K o/i ))) 5 ) 
g 4 = (l/Cl+expCcT^l.O+^Na^-^K^))) 5 ) 

g.K =A lK eX P(A K o/i) 

where 

a K =l.0,a Na = 1.0,4 = - 75 > 5 i = - 93 ^, = 2 - 6 ^ 2 = 7 - 41 > 
cr 2 = 2.0,// 2 = 2.6,<r 3 =35.7, // 3 = 1.94, 2 3 =24.3, er 4 =.88, 
// 4 = 1.4M 4 =24.6,A Na = 1.5,A 1K = 2.6,^=32.5 
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